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ABSTRACT 


A block of metal is subjected to a concentrated heat source resulting in a pool of 
molten metal surrounded by a portion of the unmelted metal. The governing system of 
equations is known as the weldpool problem. The weldpool problem is discretized using finite 
differences to discretize time derivatives and the Finite Volume Element method (FVE) to 
discretize spatial derivatives. Multigrid methods, known to be effective on uniform grids, 
make use of overlapping uniform grids of different scales to better approximate solutions 
to the weldpool problem. However, the solid-liquid interface requires extremely fine grid 
resolution and accuracy to resolve the physical behavior of the weldpool problem at this 
interface. Being too costly to apply a global fine domain, a non-uniform domain is developed 
to utilize finer resolution along the interface while still maintaining a coarser resolution on 
the rest of the domain. The fast adaptive composite grid method (FAC) is introduced, 
incorporating the concepts of multigrid to solve the weldpool problem on this non-uniform 
discrete domain. FVE is then applied to the conduction equation in the solid and the 
convection-diffusion equation in the liquid metal to develop stencil equations for use in 


FAC. 
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I. INTRODUCTION 


Several practical materials processes, such as welding and float-zone purifica- 
tion, involve a pool of molten metal with a free surface, with strong temperature gradients 
along the surface. In certain instances, thermocapillary forces are the primary drivers of 
convection in the liquid metal. In cases where other forces are stronger overall, the thermo-. 
capillary forces may still be dominant at the edges of the molten pool. Strong thermocapillary 
convection can lead to localized intense heat transfer and high velocities where the liquid 
free surface contacts the solid. This localized heat transfer can modify the shape of the 
solid-liquid interface that bounds the molten pool [Ref. 1]. 

Solving problems such as this can be challenging. One question that arises immedi- 
ately is how to discretize a continuum problem of this type. It is then necessary to construct 
a numerical process that will solve the resulting discrete system, which also poses difficul- 
ties. The purpose of this work is twofold. First, the weldpool continuum problem must be 
discretized and some numerical process developed to solve the discrete system. Methods for 
discretization will be discussed, and then incorporated into a method for finding a solution 
to the weldpool problem. Next, tools will be developed that result from this discussion. 

Several methods for discretizing and solving the system will be discussed. One method 
is the finite volume element (FVE) method, which has been shown to be an effective instru- 
ment in developing discretizations, particularly for problems that require the enforcement 
of conservation laws. Once discretized, the system of equations must be solved. Another 


technique, multigrid, has enjoyed remarkable success in solving certain types of linear and 





non-linear partial differential equations (PDEs) in a uniformly discretized domain. 

Another aspect of the problem results from examination of the solid-liquid interface. 
Near this boundary, the full souduetion-convection solution of the weldpool problem will 
require extremely fine length scales to resolve the physical behavior of the system. Multigrid 
is an efficient solver on a grid of uniform length scale, but it is not effective in cases that 
require fine grid resolution and accuracy on only a portion of the global domain. The fast 
adaptive composite grid method (FAC) is a natural extension of the multigrid methods 
applied to problems discretized on domains of varying length scales. Again, FVE serves 
admirably as the tool for spatial discretization. 

To better explain how the above techniques can be applied to the weldpool prob- 
lem, underlying ideas are first developed by applying the techniques to simple model prob- 
lems. For instance, to illustrate the method of finite differences, the time derivatives of 
the one-dimensional diffusion equation are discretized. To illustrate the finite volume ele- 
ment approach, the potential flow problem in two-dimensional space is discretized. Having 
developed an understanding of FVE discretization, the system of equations governing a sim- 
plified weldpool problem are discretized, taking into account the axisymmetric nature of the 
problem. 

FAC is then introduced as a method to deal with the need for extremely fine length 
scales along the boundary of the liquid-solid interface. FAC is an extension of multigrid 
concepts. Thus, multigrid is first developed for two model problems, a two-point boundary 


value problem describing the steady-state temperature of a long uniform rod, and the second- 


order PDE known as Poisson’s equation. With an understanding of multigrid in hand, FAC 








is developed by building on concepts from multigrid applied to a domain of varying length 





scale. Again, a simplified model problem is used to illustrate these techniques. In this 

case, the model problem is the potential fot problem in two-dimensional space. Finally, 

discretization of the simplified weldpool problem is developed for a domain with multiple 
length scales. 

The discretized equations developed are for a model of the weldpool problem, while 
the real problem is far more complicated. This work contributes to the process of solving 
the general weldpool problem by providing numerical tools for a simplified version of the 
problem. Also, complications and difficulties in the general problem may come to light by 


running the FAC discretization process on the simplified problem. 











IT. PROBLEM STATEMENT 





A semi-infinite block of metal is subjected to a concentrated heat source ap- 


plied to the horizontal surface and has an inviscid nonconducting gas above the surface. 


The result is a pool of molten metal with a free surface surrounded by the solid portion 


of the unmelted block (Figure 1). The following development of the weldpool problem is a 


recapitulation of Canright [Ref. 1]. 


The surface tension of the liquid is assumed strong enough to keep the free surface flat, 


but with surface tension variations resulting from a linear dependence on temperature, 7’. 


The resulting thermal and flow fields are assumed to be axisymmetric and steady. The time- 


dependent equations are given to facilitate a numerical approach using time-like iterations 


to reach a steady state solution [Ref. 1]. The governing system of equations is 


solid = KV°T 


oF 
Ot 
liquid : a. +u-VT=«V°T 
| 1 
OU ea Vase vos 
Ot p 
V 


solid surface z=0 : as = 0) 
Oz 
liquid surface z=0 : ped = —q(r) 
Oz 
v= 0 
Ou _ Oo 
Moz ar 


(11.1) 
(11.2) 
(11.3) 


(II.4) 





axis r = 0 a0 8 
u=0 
Ov 
aay 
Tr 


This system is governed by conservation of energy in the solid and by conservation of energy, 
momentum, and mass in the pool. In the system, « is thermal diffusivity, t is time, p is 
density, p is pressure, v is kinematic viscosity, u is velocity with components u and v in 
the r and z directions (cylindrical coordinates), k is thermal conductivity, g(r) is imposed 
surface heat flux (large at r = 0 and falling off to zero at some small value of r such that 
So° a(r)2nrdr = Q), p is viscosity, and y is the negative of the derivative of the surface 


tension with respect to temperature (y is assumed to be constant and positive). 


Q 


insulated r 





solid 


Figure 1. A block of pure metal is subjected to a concentrated surface heating, Q, that results 
an a molten pool. 








The equations can be nondimensionalized. The main dimensionless parameters are 
the Marangoni number Ma and the Reynolds number Re. It is also convenient to eliminate 


the pressure by adopting a stream function/vorticity formulation for the flow [Ref. 1]: 


Re (Fv x (uxs)] = -VxVxw 
gC = vx vx (4) 
u=V x (=) ¢ = DOE ee 
r r Oz r Or 


where W is the axisymmetric stream function and w is the vorticity vector (having only one 


component, in the 6 direction). The flow boundary conditions are 


liquid surface z=0 : Y=0 
se OE 
Or 
axisr=0 : V=0 
w = 0 
o. ; OV ow 
liquid/solid boundary r = f(z, t) Y= rar 


All of this results in the following equations: 





solid : ad = Ma !W°T (II.5) 

— oT ~1—02 
liquid : > + V -(uT) = Ma V°T (11.6) 
SV x (uxw) = —ReIV x V xw (II.7) 

W\ . 
w=VxVx (=) 6 (11.8) 
W\ « lov. 10Wv. 

where u=V x (=) d= ee = = aes (11.9) 





with the conditions at the boundaries and at the solid-liquid boundary given as 


solid surface z = 0 


liquid surface z = 0 


axis rT = Q 


Refer to Canright [Ref. 1] for a more in-depth development of the weldpool problem. 


This system of equations and systems like it present interesting challenges for con- 


structing numerical solutions. The goal is to transform this continuum problem into a 


discrete one and to develop a numerical process to solve it. The Finite Volume Element 


Method (F VE) has proven to be a useful tool in discretizing systems involving the enforce- 


ment of conservation laws. Multi-level techniques, such as multigrid, have been effective in 


streamlining the solution of partial differential equations (PDE’s). The Fast Adaptive Com- 


posite method (FAC) is designed to efficiently discretize and find solutions for PDE’s. FAC 


methods are characterized by their use of a composite grid, the union of various uniform 


grids. FAC methods require discretization of the PDE on both the individual levels and the 


composite grid. F'VE is utilized for this purpose. 














Ill. DISCRETIZATION 


Numerically solving the conduction problem 


= . Iif.1 
Ai KV*°T (III.1) 
surface z=0 : co = —d(r) 
axis r=0Q : a = 0 
far boundaries r,z=1l : T=0 


requires discretization of both the space and time derivatives, where the value of the unknown 
function T is determined from a set of equations which approximates equation (III.1). In 
the following, finite differences are used to discretize the time derivatives, while the Finite 


Volume Element (FVE) method is used to discretize the spatial derivatives. 


A. FINITE DIFFERENCES AND CRANK-NICOLSON 
Discretizing in both time and space can be complicated, and it may be worthwhile to 
develop the method on a simple one-dimensional equation. To illustrate the discretization 


of the time derivatives, consider the one-dimensional diffusion equation, 


= ae (III.2) 


Assuming the function T and its derivatives are single-valued, finite and continuous functions 


of x and t, then Taylor’s theorem can be used to approximate each term of (III.2). Holding 


t constant (i.e., treating J as a function of x only) and applying Taylor’s theorem gives 





OF. dn gO 1 .. Ly fOr 
a se ET A, eg ene ag a, eS ITI.3 
T(z + Az) T(z) + Ara + 5A age t eo a3 t (II1.3) 
and 
2 
is ais eA Neda: (III.4) 


=a — =a nme _ ee ~— 
a—Ag) ]T (x) Ar + Aros Ae ~ GO? ps 


Summing these expressions gives 


2 
T(z + Az) + T(x — Az) = 2T (xr) + Arte + O(Ax"*), 


where O(Ax*) denotes the terms containing fourth and higher powers of Ax. These terms 


are assumed to be negligible compared to the lower powers of Az. This gives 
OT 1 
ae Ae [T(z + Az) —2T(z) + T(z — Az). (III.5) 
Let 7; represent the temperature at time nAt and spatial location kAz. Using this notation 


Tp = T(kAz,nAt). Then, in two dimensions, equation (III.5) is written as 


OT es 14 a 21, + Tey 
Ar? ~~ Ax? 


Other derivatives can be approximated in similar fashion. For example, subtracting 


equation (III.4) from (III.3) and neglecting the terms of order Az? and higher gives 


(III.6) 


a A T(e + Az) — T(x - As). 


ae 


Equation (III.6) is called a central difference approximation. 


For the time derivative (holding the spatial term constant), OT /Ot may be approxi- 


mated in several ways. The forward difference is given as 


OT Tpt'-TP 
Ot i 


10 











The backward difference is 


The central difference is derived in a manner similar to that leading to equation (III.6). This 


results in 


OT dae a dae 
fy yA. rn 
Using the forward finite difference approximation, the one-dimensional diffusion equa- 


tion (III.2) can be approximated by 


Tet Te _ Tha — 2Te + Tes | 
oor ah , (111.7) 


Letting r = At/Az’, this equation can be written as 
Tet! = T™ — (1—2r)TR +r Tey. (III.8) 


The unknown 77+! can now be computed explicitly at the (n + 1)st time step using the 
known values of the nth time step. 

While this explicit relation provides a simple means to compute unknown values, 
the process is only stable for 0 < At/(Az)? < 1/2 [Ref. 2]. By imposing a stability limit 
on At/(Azx)?, explicit finite difference methods for the diffusion equation can give useful 
approximations. However, the requirement that At/(Az)* < 1/2 places a severe restriction 
on the interval size in the t direction resulting in an increase in computation. A better 
method is needed. 

An implicit formula is one in which two or more unknown values in the (n + 1)st 


time level are specified in terms of known values in the nth time level by a single equation. 


11 





If there are M unknown values in the (n + 1)st time level, the formula must be applied M 
times across the length of the time level. The resulting system of M simultaneous equations 
specifies the M values implicitly. A method from Crank and Nicolson [Ref. 2] combines a 
central difference approximation to the time derivative at the half time-step (n + 4) with 


the average of central difference approximations to cae at times nAt and (n+ 1)At to 


aT n+i OT n+4 
(J | = (Fa) | (II1.9) 


k k 


approximate the equation 


This equation is then approximated by 


Tet Tp 1 (Tee OTP + TEN | TR — P+ TR (111.10) 
At 2 Az? Ax? 
This reduces to 
—rT py + (24+ ar)Iet! — rT = rT, + (2-2r)TP + rTP, (III.11) 


where, again, r = At/(Az)*. There are now three unknowns at the (n + 1)st time step 
written in terms of three known values at the nth time step. 

Although Crank and Nicolson’s method is stable for all positive values of r, large 
values of r, e.g., 40, have been shown to introduce oscillations in the solution (Ref. 3]. To 
combat this problem, a more general finite difference approximation is implemented, using 
a weighted average of the finite difference approximations to oF at the nth and (n+ 1)st 


time steps. This weighted average is given by 


Tyt!-— Tp Tei — Tet + Te Te, — 278 +TR. 
WATE | DAT ed 1— ae on eed III. 
ie a ( Aa +(1-a) Ag ; (III.12) 


where 0 <a <1. Notice that a = 0 gives the explicit scheme, a = 1 gives a fully implicit 


backward difference method, and a = 1/2 gives the Crank-Nicolson method. This weighted 
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average is stable and convergent for all r if 1 < a < 1, but is stable and convergent for 
0<a< i only when r < 1/(2 — 4a) [Ref. 3]. 
Following the method outlined above, a general finite difference approximation for 


the time derivative of the two-dimensional conduction problem (III.1) 


= = WT (111.13) 
is given by 
m+l _oyn 
ss = akV2T"™*! + (1 — a)KV°T”. (III.14) 


In this approximation, 7” is the temperature at the nth time step. The current time step 
is designated by n and the next time step by (n + 1). Combining terms for the current and 


subsequent time steps yields the semi-discrete relationship, 


22 a 4 m+1il _ i Ee 2 
(x akV )r = (=; +(1—a)kV i tee (III.15) 


The above conditions on a also apply here. 


B. FINITE VOLUME ELEMENT METHODS 

Having dealt with the time discretization, it is now necessary to address the dis- 
cretization of the spatial derivatives. This development of discretization for spatial deriva- 
tives closely follows McCormick [Ref. 4]. The spatial portion of the problem is discretized 
using the Finite Volume Element method (FVE). The classical finite volume method (FV) 
is commonly used as a discretization method for computational fluid dynamics applications 
because of its ability to be faithful to the physics in general and conservation in particular. 


However, use of F'V requires a scheme for approximating certain fluxes, which is often done in 


13 








an effective but rather ad hoc and restrictive way. Thus, the Finite Volume Element method 
was developed in an attempt to use the flexibility of the finite element representation of the 
oe functions to create a more systematic FV methodology [Ref. 4]. The central idea 
is to approximate the discrete fluxes needed in FV by replacing the unknown PDE solution 


by a finite element approximation. 





Figure 2. Control volume V (dashed lines) and surface S. 


As an example, consider the potential flow problem in Euclidean two-dimensional 
space: 
V-(pVv) = n in (III.16) 
w = yo on dOQn UNE (Dirichlet) 
(pVwv)-nrm = Y on OQwUOQs (Neumann). 


The domain for this problem is the unit square 2 = [0,1] x [0,1] € R*. In this domain, 


OQn and OOg are the northern and eastern Dirichlet boundaries of the domain 2, and 005 
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and OQy are the southern and western Neumann boundaries. Beginning with a “control 
volume” V in the domain 2 with surface S (Figure 2), the potential equation (III.16) is 


integrated over the volume V to get 


I V-(pVw)dV = [ ndV. (III.17) 


Using the Gauss Divergence Theorem, the left-hand side volume integral is transformed into 





Figure 3. Volume partition (dashed lines) of Q. 


a surface integral, 
[ (oVu) -AdS = I ndV. (III.18) 


Because this problem is two-dimensional, “volume” refers to “area” and “surface” refers to 
“boundary”. In other words, volume integrals in two dimensions are area integrals, and the 
surface integrals are line integrals taken over the bounding curves. 

Each side of equation (III.18) represents a flow rate in mass per unit time, and 


(pVw) «7 represents the flux density across the surface 5. The integral condition (III.18) 


15 





can be interpreted as a conservation law for the volume V. Conservation can be enforced 
over the domain by requiring (III.18) to hold over the individual volumes. The potential 
equation (III.16) can then be discretized by choosing a finite set of volumes that partitions Q 


into volumes (Figure 3) with the integral condition in (III.18) being imposed on each volume. 


Figure 4. Finite element partition of the domain QQ. 


It is now necessary to select the functions to be used to approximate w. It is well 
known that piecewise linear functions are often effective at approximating functions in finite 
element methods. Thus, VW can be approximated by continuous, piecewise linear functions. 
Then, a basis — a set of functions whose linear combinations determine the desired space of 
functions — is selected. The directions of interest are x and y with a uniform grid of step size 
— + applied in both directions. Each grid square is divided into two triangular elements 
by a diagonal oriented in the direction of increasing x and y (Figure 4), and “hat” functions 


are constructed using these elements. The FVE method replaces the exact solution * with 
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a finite element approximation yw expressed as linear combinations of the basis functions. 
At each grid point, the value of the approximation equals the weight for the associated hat 


function. Thus, ~ is expressed in terms of a basis: 


N 
p= \~ Ui j Pig (Zs Y)s (III.19) 


ij=l 

where double subscripts, 7,7, varying from 0 to N, are used to describe the location of the 

grid point (z;,y;). Then, u;,; is the value of y at the (2,7)th grid point, (z;,y;), and ¢,, is 

the “hat” function, a continuous piecewise linear function associated with the (i, 7)th grid 
point. 

To simplify notation, the sampled values of the unknown w, currently written as 

an (N + 1) x (N +1) two-dimensional array, can be considered a one-dimensional array. 

To accomplish this, the elements 7; are labeled in lexicographic order. For example, the 


elements in the two-dimensional case are written as 


Yn Una ++ YNN 


P10 W1,1 = Wi,N 
Wo,0 Wo,1 anaes Wo,N 


These elements can be written as 


T 
(wWo.0, Wo,1) yes Yon: Y1,0) Y115 cy WIN; cee » YN; Wn, cer y wn.n) , 


This array is normally considered to be a column vector. Using this notation, ~ is expressed 


in terms of a relabeled basis: 


(N+1)? 


y= a uoy(Z, y). (IIT.20) 
t=1 
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Figure 5. Hat function associated with the (1, 7)th grid point. 


Substituting (III.20) into (III.18) and integrating over each volume V, (following the 


same ordering scheme) for k = 1,2,...(N +1)’, the left side of equation (III.18) becomes, for 


each k, 
(N+1)? 
[eve dSe= ou fh oe Vor- reds. (111.21) 
Sk l=1 St k 
This gives the system Lu = f, where 
f= | ndVz (III.22) 
Vi 


(except on the boundary) and L is the (NV +1) x (N +1) matrix with entries 
Det = [ (p, Vo) ; Npd Sy. (III.23) 
k 


To correctly implement the boundary conditions, the system must be modified. Since 


both Neumann and Dirichlet conditions will be of interest, both cases will be examined. 
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Neumann conditions are imposed indirectly on ~ by substituting the flux value ~, from the 
boundary conditions into the appropriate term in equation (III.18). The Neumann conditions 


in equation (III.16) give 


[ OV Ads = [ ys. BS 4 [ syseOVe) fds = I ndV, (II1.24) 


where Sw and Sg represent the southern and western bounding curves of the volume where 
Neumann conditions are imposed. 

By contrast, Dirichlet conditions in equation (III.16) are imposed directly on w. 
Thus, ~ takes on the value of yp at each Dirichlet grid point, or rather each grid point 
on OQw UONeg. There are potentially fewer 7 unknowns than equations since there are now 
volumes that do not have unknowns. This is easily avoided by expanding the volumes at the 


points neighboring the Dirichlet boundary (Figure 6). 


Figure 6. Modified volume partition of Q to accomodate Dirichlet points (e). 


In addition, the intersection of OQy and ONw and the intersection of Os and OQz 
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are treated as Dirichlet grid points. These points could also have been treated as Neumann 
grid points. For full Neumann problems, the original partition of Figure 3 is used. 

It is now time to evaluate the integrals in (III.22) and (III.23). For the integral in 
(III.23), the rule 


[ (eve) -AdS = p(Pr) A Vo: nds (III.25) 


is used on each bounding segment, T, of the surface. That is, T is one of Sn, Sg, Ss or Sw, 
and Pr is the point of intersection of T and the grid lines passing through Ny, the grid point 


at the center of the volume. For Neumann boundary segments Ss and Sw, the rule is 


[vids = vi (Mr)ITI, (111.26) 


where Mr is the midpoint of T and |7| is the “surface area,” or length, of T. 

Consider the discretization on a uniform mesh with grid size h = 7 in both coor- 
dinates. Returning to the double subscript notation, 1 and j7 vary from 0 to N — 1, where 
0 corresponds to the Neumann boundary grid point, or node, along OQy and OQs5, and 
N —1 corresponds to the nodes next to the Dirichlet boundary, Oy and ONg. For a general 


interior node (i > 0,7 < m-— 1), consider one segment of the boundary, say the NW side 


(Figure 7). To evaluate the integral 
o(Pr) : Vw -ndS, (I11.27) 
NW 


first notice that, with n = ( ; ) 


Tae = oe 


; (III.28) 


since 7~ is piecewise linear. Then, 
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i,jel 141,341 





i-1,j-1 i,j-1 


Figure 7. Volume of general interior grid point (i, j). 


a0. Pega — Vig _ A (Vigir — Vig) _ Vigo — Yi 
i= Vy fads = itt — ‘ dS = 5 (| Se  (L29) 


This process is repeated on all segments of the volume boundary (8 calculations in all). The 
NW and NE calculations can be combined into a single North calculation. The same occurs 


for the South, East and West sides: 


I Vi-adS = Viger — Vig (III.30) 
[ Vb-adS = bia Vij (111.31) 
[ Vi-AadS = bing — Vig (III.32) 
[Ve nas 2 peasy (11.33) 


Then, for the general interior grid point, or node, 
p(Pr) [ Vy nds = o(Pw)i-1g +e( Pe) Vit1g + e(Pn) Vig te(Ps)vig-1— >, Vig, (II1.34) 


where > = p(Py) + p(Ps) + p(Pe) + (Pw). 
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The quadrature rule is used on the integral in (IIJ.22). This results in 
/ ndV = nig|Vj\, (III.35) 
Vj 


where 7; is the value of (x,y) associated with V; and |V;| = fy, dV is the “volume,” or 
area, of V;. For a general interior node, |V;| = h?. Thus, for a general interior node, the 
stencil equation is 
Pij+ 
P34; D Pisd yj Wij = h’n(ih, jh), (111.36) 
Pig-t 


where »& signifies the sum of the outer stencil entries and p;,; = p(th, jh). This notation for 


p(ih, 7h) will also be used in subsequent stencils. In the stencil notation, 

ABC 

DE FY 45, 

GH I 
stencil entries (A, B, ...) indicate the value to be applied, with the position of the entry 
indicating to which grid point the value is to be applied. Specifically, the center of the 


stencil is the current grid point, with the surrounding stencil entries corresponding to the 


neighboring grid points in those directions on the grid. Thus, the above stencil has the value 
Ai-1j41 + Bui jar + Cbinr ga + DWi-1j + Evig + Puieig +.... 


For a side Neumann node, stencil development follows a similar path. Along the 
boundary, recall that the Neumann condition gives Vy-:n = yy. Integrating on the Southern 


boundary (Figure 8) and using the midpoint rule, 


[ Vila, y)dS = heh (2%, 0). (11.37) 
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For the smaller volume, 


h? 
: ndV = nij\Vi| = ni (=) (III.38) 


For a side Neumann boundary grid point (j = 0 on southern boundary), the stencil is 


he 2. 
2Pi-4o = > SPixd o | Yio = meh 0) — hy, (ih, 0). (III.39) 


2? 


For the corner Neumann node (i = j = 0), the stencil is developed in the same manner as 


amqnr nmr wr er re — a a a: on 





i-1,0 1,0 i+1,0 


Figure 8. Side Neumann grid point on the southern boundary (j = 0). 


the side Neumann nodes. The stencil is 


1 
2Po,$ h? h h h 


The stencils for the neighboring Dirichlet grid points (i = N — 1 or j = N — 1) are 
the same as at interior grid points, except that the Dirichlet points are eliminated and the 
stencil entry reaching to a Dirichlet node value is removed and placed on the right-hand 


side as a coefficient of the boundary data. There are two possible stencils for a neighboring 
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Dirichlet node: the reduced volume and the extended volume. For the neighboring Dirichlet 


nodes of the Eastern boundary, OQ, (t= N-1,0 < 7 < N—1), the reduced volume stencil 
1S | 
PN-1,j+2 


Py-ag DL 


Wn-1,; = hén(t i h, jh) = Py—2 soll, 7h), 
PN-1,j-4 


(III.41) 
where » is the sum of the outer stencil entries before the boundary terms are removed. For 


this case, © = p(1 — 4 


33h) + p(l— By jh) + p(1 —h, (J — 4)h) + o(1 — A, (7 + 4)h). 


The stencil for the same node in an extended volume is 


1 
2PN-1,j+4 
3h? 1 is 
PN-3j = >», WnN-1j3 = > (1 —h, jh) + 5Pn-3,5¥o(1, gh) 
Pn-1,5-4 


—pw,jvo(1, (7 + 1)h) - =px.so(l, (j — 1)h) 


where & is, again, the sum of the entries of the stencil without the boundary terms removed. 


The extended volume stencil for the neighboring corner Dirichlet point (i = N —1 and 
j=N-l1)is 


9h? 
3PN-3,N-1 = > 


i 
WN-1,N-1 ras a Sy A h) ~ 5 Pn -4,n—1¥0(1, l= h) 


3PN-1,N-3 
1 
~Pw-4,-3¥o(1, 1 — 2h) — Spy_1n—4¥o(1 — hy 1) 
—Pv—-3,n-i¥o(1 = Zh. 1). 


C. DISCRETIZATION OF THE WELDPOOL PROBLEM 


The technique used on the potential problem (III.16) may be used to discretize the 


weldpool problem. To begin, the FVE stencils for the conduction problem (II.5) are devel- 
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Figure 9. Toroidal volume for the conduction problem. 


oped. For the axisymmetric conduction equation, a control volume for the lth grid point is a 
toroidal prism. This volume, V;, is generated by taking a two-dimensional volume for the /th 
grid point in the (r, z) plane and sweeping the volume through 360° in the ¢ direction (Figure 


9). Previously, a general finite difference semi-discrete approximation for the equation 


a = Ma !V’T 


was found to be 


Trt — 7? 


oe aMa!V?T""! + (1-a)Ma?*V*T". 


Combining terms for the current and subsequent time steps yielded the semi-discrete rela- 


tionship, 


ei -1~02 n+l __ ji _ —ly2 n 
ts aMa v’) 7 =(4+0 a)Ma“V*)|T". (III.42) 


To complete the discretization, discrete equations for the spatial derivatives are developed. 


‘ The resulting equations are assembled into an algebraic system. Begin by integrating (III.42) 
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over all of the control volumes to yield 
1 —~1~02 ) nti / 1 Vf .-102 
eee ald = am ma "dV. III.43 
La aMa vi) dv = | (—+(1-a)Ma'V?)T (III.43) 
Here, V represents the partition of the domain 2 into the smaller control volumes, V;, or 


V =UYV,. Applying the Gauss Divergence Theorem gives 
=: | T™ dV —oMa"? / VT™ ads = / T"dV +(1-a)Ma7 | VT" nd, (II1.44) 
At Jv S At Jv s 


where 7 is the outward normal. 

As in the mode! problem, 7 can be approximated by continuous, piecewise linear 
functions. Using cylindrical coordinates, a basis is chosen for the space in which these 
functions are found. Because of cylindrical symmetry, the discretization can be examined 
on the half plane 6 = 0. Hence, T is treated as though it is a function defined on a two- 
dimensional domain in the (r,z) plane. In this two-dimensional domain, “hat” functions , 


¢i3(7, z), are chosen as the basis functions, and the unknown function T is approximated as 


N 
L(h2) = s~ 1550; 4 (72), (III.45) 


2,7=0 


where the uniform grid step size is h = 1/N. 

A change of notation is again used. The sampled values of the unknown T, now 
written as an (NV +1) x (N +1) two-dimensional array, can be considered a one-dimensional 
array. Io accomplish this, the elements T;; are labeled in lexicographic order from | = 1 to 


1 =(N+1)*. For example, the elements in the two-dimensional case are written as 


Ino Ini «+: Tn 
Zio fia +++ Tin 
Too Loi ‘*: Town 
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These elements can be written lexicographically in a column vector as 


fe 
(Too, To,1; ve ,10,N; Ti0; Ti 1; vey TiN; a Two, 2,1; nee , Inn) 


The values on the grid 92 are given now as 
Ti = T (ri; 21). 


Substituting equation (III.45) into equation (III.44) gives 


(N+1)? 1 
n+l 
7, Fe i 


eptdv — aMa” [ Vent! ads| (111.46) 
(N+1)? | 


= » TT x [sav +a-o) )Ma- | var -ads|. 


The volume V of the domain is partitioned into 1 = (N + 1)* volumes (actually, there are 


fewer volumes because of boundary conditions) where 
V=WN and Vif )v=0 for LFk. 


Substituting the partitioned volumes into equation (III.46) results in 


(N+1)? 


Trt 
2, 


: (N+1)? (N+1)? 


— ae grav, —aMa"! x [ V¢rt!. ads, (II1.47) 





(N+1)? (N+1)? (N+1)? 
= 7 x y [ stay. + (l-a)Ma" YO / Vor -AdS;,| . 
i=1 | k=l “9k 


Notice here that the interior surface integrals cancel out. The eastern edge of one interior 
volume coincides with the western edge of the neighboring volume. This edge is integrated 
over twice, but in opposing directions, resulting in the cancellation. 

This set of equations can be written as an operator equation, L[T"*'] = f(T"), or as 


a matrix equation, M,[7"*}] = M.[7"]. The operator L and the elements of the matrix M, 
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are given as 


1 -lw2 
a se 111.48 
L : ( x - Mav dV, (III.48) 
Mig =e onlay, —aMa™! / Vor! . adSp, (III.49) 
At Vz Sr 


while f(7™) and the elements of Mz are given as 


1 
f(T") = I [+0 - a)Ma'v* T"dV (III.50) 
1 n _ no 
Mu = = z bpd, + (1 —a)Ma"! [ Vol Ads, (111.51) 


Any function whose coefficients satisfy the resulting set of discrete equations also 
satisfies the conservation law over any volume made up of the union of control volumes. 
However, complications arise when considering the boundaries. The Neumann conditions 
on the western and southern boundaries are incorporated indirectly into T by substituting 
the normal derivative value into the appropriate term in equation (III.44). The Dirichlet 
conditions on the northern and eastern boundaries are imposed directly on T. The control 
volumes are extended so that the control volumes normally associated with the boundary 
points are taken into account by the extended volume of the neighboring points (Figure 
10). The integrals over the basis functions are calculated by representing the basis functions 
as collections of triangular planes. These triangular planes are assembled to form the hat 
function (Figure 5). An interior grid point has six triangular planes joined in this manner, 
each triangular plane corresponding to the six triangles surrounding the grid point (Figure 


7). Integration results in the following stencils for a general interior point on a uniform grid: 


2 o2—0€ 16+ 5¢€ 


h 
| (x: Toate dV = — | 32-ile 224 32+11e | 7. 
LA Pq 384 ig be 32 i - 
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Figure 10. A cross section in the (r,z) plane of a finite volume element grid with the domain 
partitioned into control volumes (dashed lines). Extended control volumes are used for grid 
points neighboring Dirichlet points. The triangular finite elements (solid diagonal lines ) form 
the piecewise linear “hat” function centered at each node. 


| 2 
; (So ToaVane 8) dS = 5] 2-« —8 2+€ 


Pq 2 





where € = & (r is the radial coordinate for the central point of the volume). The stencils 
use the two-dimensional labeling. The summations signify summing over all “hat” functions 
for each volume. Also, the volume and surface integrals are in cylindrical coordinates. This 


means that the integrals are as follows: 


[av = [ rdrdodz = 2n [ rdrdz 


[ ds = [ rdOdl = 2n $ rd. 


In each case, the 27 resulting from the integration in the @ direction has been factored out. 


The discretization of the convection-diffusion equation, 


a +V+(uT) = Ma'V’T, 
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has the same form as the conduction equation for the solid except for the term 


[ Vv: (uT)dV 


Using the Gauss Divergence Theorem, this portion of the equation is also discretized. By 


cylindrical symmetry, 
/ V -(uT)dV = [ (uP) fidA = 2 § (uP) - Aardl, 
V 


where the 27, as above, is factored out. Remembering that 


_ tl OW . 1ow 5 
 r Oz r Or 
produces the result 
g(ul) -nardl = [Tar - f _ <—Tdr r+ fF oe Te — f Tae. 


Substituting the piecewise linear “hat” functions for the unknowns WV and T into the above 
integral results in discrete equations. Evaluating and summing the integrals and using stencil 


notation gives 


[.v-(arav = 


; Y =| “i 
; -1 |w fi 1 | <i Ww | T;, 
—] =| y 
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Discretization of the vorticity (equation II.7) follows the same scheme. Integrating 


the vorticity equation over a control volume results in 
=f i. dV -| Vx (ux a)dV = -Re™ | Vx V xadV. (III.52) 
Ot Jv V V 


The integration over the volume equates to an integral over the cross sectional area with 6p 


| acting as the normal vector. Applying Stokes’ theorem to (III.52) results in 
2 / / wdrdz — $ i. (ux d)dl = —Re™ $ £.(V x d)dl. (111.53) 
Ots JA C C 


Again, a uniform grid with step size h in both the r and z directions is applied. Each square 
of the grid is divided into two triangular elements by a diagonal in the direction of increasing 
r and z (Figure 10) and linear interpolation is used over each triangular element. The control 
volume cross sections are square with side length h and centered on a grid point, except for 
the control volumes surrounding side and corner Neumann grid points and side and corner 
Dirichlet boundary points (again, see Figure 10). 

The area integrals are taken over six triangular regions, while the line integrals are 


taken over four line segments, each line segment containing halves in two different elements. 


Recall that 


10ov¥. lowv. 
u=—~7t + a4, (III.54) 


and note that the vorticity vector, ©, has only one component, the component in the 6 


direction. Using the expression for u from above results in 


. wOoV wodV. 
ux @= —-——f -— —-—z. 


r Or r Oz 


The integrals in (III.53) then reduce to 
ov w OV w OV w OV w 
af fjwdrds + f Star [ ocars [Zo Sde— [ode 
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= Bee (fear - [ Sear- f ed ae + [ eae) 


After substituting the piecewise linear element representation of the unknowns into 
the above integrals, the discrete equations for an interior grid point on a uniform grid are 


given in stencil form as 


2 
set] —l de 
—2 1l- set Ww 1 |W 
1 2* 1+ —1+3er €~ —l~ 
aa set —l|—et Y l+eT -eT+e lLl-e |V —l+e -fe7 | U Iw 
—i+ —«rt —l— He 2*1-— 


where € =h/r, e* =e/(1—€/2), € =e/(1+6e/2), 17 =1/(1—/2), and 17 = 1/(1+€/2). 


The last equation to be considered for discretization is the stream function (II.8), 
w=Vxu4u, (111.55) 


where 


Integrating over the control volumes and applying Stokes Theorem to the right side of equa- 


tion (III.55) gives 


/ / wdrdz = g f- udl, (11.56) 
A C 
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where # is, again, the unit vector tangent to the curve C. With u as above (equation III.54), 


these integrals reduce to 


| [ wdrdz = eee ear + | oe nde - f ce (111.57) 
A Ww E 


IN Oz r 5 Ozr Orr Orr 


Again substituting the piecewise linear element representation of the unknowns into 
the above integrals gives the discrete equations for an interior grid point on a uniform grid. 
In stencil form, 
dnr{. 2 } | = ; 

——| 2 142 ]w = -1t (2+1°7+4+1*) -1 Wi5, _ (TII.58) 
1 —1 
where the values for 1~ and 1* are as above. 

The weldpool problem is now discretized on a uniform grid. In the next section, the 

discussion will center on what a composite grid is, why a composite grid is needed, and how 


to resolve the weldpool problem more efficiently. 
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IV. FAST ADAPTIVE COMPOSITE METHODS 


This chapter will introduce the concepts of fast adaptive composite (FAC) 
grid methods. FAC methods are used to discretize and find solutions for partial differential 
equations (PDEs) on a composite grid. A composite grid is formed by the union of a nested 
sequence of uniform grids of varying size and scale. A patch, or rectangular uniform grid, 
can be aligned with coarser grids and can be épecified by its relative location on the coarser 
grid, its dimensions, and its mesh size (Figure 11). The PDE is discretized and solved on 
this composite grid. However, the actual computations take place on the uniform subgrids. 
In this manner, FAC maintains the advantages of uniform grid discretizations ate allowing 
effective adaptation to the local phenomena. FAC methods require discretization of the PDE 
on both the individual levels and the composite grid. The approach used here, the finite 
volume element (FVE) method, combines the finite volume and finite element methods (Ref. 
4]. 

Many physical models encountered provide excellent opportunities to employ 
FAC. Often, physical models require higher accuracy only in limited regions of a domain. 
The physical model could be resolved on a global uniform grid with the smallest required 
mesh size, but computational power is rarely unlimited, and some method for locally adapt- 
ing resolution and accuracy is needed. It can also be an advantage to introduce as little 
nonuniformity into the method as possible. This is accomplished using composite grids. For 
an example of a coarse grid, consider the two-level grid in Figure 12. The coarse grid is 


defined on a domain {2,, where 2h indicates the grid spacing. There are several distinct 
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Figure 11. A composite grid with its fine grid patch. 


types of gridpoints, such as regular grid, coarse grid interface, fine grid interface, and slave 
points. The fine grid points can be treated as regular interior points of the refinement patch, 
while coarse grid interface points and slave interface points can be taken as boundary points 
for the fine grid patch. Since the functions on all of the grids are piecewise linear, the slave 
interface points can be calculated as a linear interpolation of neighboring coarse grid in- 
terface points. Multilevel methods can then be implemented. These methods make use of © 
fully overlapping grids of different scales and correction of the coarse grid equations by the 
composite grid residuals to approximate the PDE solution. This is accomplished by solving 
the system on the fine grid level and then correcting the coarse grid equations with the fine 


grid approximations. 
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& fine grid interface ‘a slave pointe C) coarse grid interface 


Figure 12. Composite Grid Interface Points. 


A. FUNDAMENTALS OF MULTIGRID 

Before entering into a development of FAC, the concepts of multigrid should be in- 
troduced. Multigrid is a basic component of FAC because of its apparent effectiveness as 
an iterative method. It is of interest to determine how multigrid behaves as a solver of 
discretizations on uniform grids, because multigrid will be used for composite grid equations 
restricted to the global and local subgrids [Ref. 4]. 

Multigrid methods were first developed to solve boundary value problems posed on 
spatial domains. These boundary value problems can be made discrete by selecting a set of 
grid points in a domain. The discrete problem that results is a system of algebraic equations 
associated with the chosen grid points. The simplicity of these boundary value problems 


provides a good introduction to multigrid. 
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1. Two Model Problems 
Consider the two-point boundary value problem describing the steady-state temper- 
ature distribution in a long uniform rod. The second-order ordinary differential equation 


1S 
=u (t) euls) = fle); 0<a< 1h. oS 0, (IV.1) 


where u(0) = u(1) = 0. The simplest method for solving this equation numerically is the 
finite difference method. The domain of the problem is partitioned into N subintervals in a 


step size of h=1/N. A grid is established with points 2; = jh (Figure 13). 





i 


h 


Figure 13. One-dimensional grid on the interval 0 < x < 1 with grid spacingh =1/N and 
grid point x; = gh forO<j <N. 


At each of the N —1 interior grid points, the differential equation (IV.1) is replaced by 
a finite difference approximation. Since the exact solution is unknown, the approximation v, 
is introduced. The approximate solution is represented by the vector v = (v, v2,-.-,UN-1), 


whose components satisfy the N — 1 linear equations of the finite difference approximation, 


Ua a UU 


pT tou = F(t), ISiSN-1, (IV.2) 
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with the boundary conditions v) = uy = 0. Representing this system of linear equations in 


matrix form gives 


2+ ah? = | Vv] fi 
—-l1 2+o0h? -1 
ip a wi Tet. (IV.3) 
; _] : 
—1 2+ ah? UN-1 fn-1 


or simply Av = f. The matrix A is tridiagonal, symmetric positive definite, sparse (the 


matrix elements are predominantly zero), and has dimension (N — 1) x (N — 1). 
For the two-dimensional case, consider Poisson’s equation, a second-order partial 


differential equation given by 
—Uge — Uyy = f(z,y), O<zr<1, O<y<l, (IV.4) 


with u = 0 along the boundary of the unit square (Figure 14). 


Y 


Xx 


Figure 14. Two-dimensional grid on the unit square. 
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As before, Poisson’s equation can be cast in a discrete form using the grid points 
(xi, y;) = (ih, jh), where h = 1/N is the width of the subintervals. Replacing the derivatives 
of Poisson’s equation (IV.4) with second-order finite differences and the exact solution with 


the approximation v;,; gives 


~Vinig + vig — Vieng | THT 20ij — Vij-l _ aa, (IV.5) 


h2 h2 
with uv; = 0 along the boundary (i=Oori=Norj=0orj=N)andl1<i<N-—1Jand 
1 <j< N-1. There are now (N — 1)? interior grid points and unknowns. By arranging 
the lines of constant i in lexicographic order, the unknowns of the ith row of the grid can 
be collected in the vector v; = (vj1, Vi2,---,Ui,n—-1)’ for 1 <i < N—1. When f (xi, y;) is 
ordered in the same manner, the system of equations (IV.5) can then be written in matrix 
form as 

A —I Vi fi 
-—. A -I 
l , : 


< eae (IV.6) 


= VN-1 fy—1 


This system is symmetric, block tridiagonal, and sparse. Each block on the diagonal 
is an (N — 1) x (N — 1) tridiagonal matrix, while each off-diagonal block is the negative of 
the (N — 1) x (N — 1) identity matrix J. 

2. The Residual Equation 

Remember that a system of linear equations, such as equations (IV.2) and (IV.5), 


can be written as 


Au =f ; (IV.7) 

















with u being the exact solution. Let v represent an approximation to the exact solution. 
The exact solution is unknown, while the approximation v, which is often computed using 
an iterative method, is known. One important measure of v as an approximation of u is 


algebraic error, given by 


e=u-Vv. (IV.8) 


However, since the exact solution is still not known, the error is also unknown. Still, one 


measure of how well v approximates u is the residual, given by 
r=f— Av. — (IV.9) 


The residual is the amount by which the approximation fails to satisfy the original problem 


in (IV.7). Rewriting the residual, (IV.9), as 
Av=f-r (IV.10) 


and subtracting this new version from the original system, equation (IV.7), produces the 


residual equation 


Ae =r. (IV.11) 


Assuming that the approximation v has been computed by some method, the residual r can 
be easily computed using r = f — Av. To improve the approximation, solve the residual 
squation (IV.11) for e and compute a new approximation using the definition of the error 
(IV.8). This is the idea of residual correction and will prove to be very important (Ref. 5]. 
The residual equation shows that the error satisfies the same set of equations as the 


desired exact solution when f is replaced by r. This suggests that performing some iterative 
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method on the original equation Av = f with an arbitrary guess v is equivalent to performing 


the same iterative method on the residual equation with the specific initial guess e = 0. 


3. Solving Systems of Linear Equations 

There are several existing methods for solving sparse systems of linear equations 
like equations (IV.2) and (IV.4). Direct methods solve the system of equations in a known 
number of arithmetic operations. Several types of errors can result from using direct methods. 
Algebraic errors in the solution result from rounding errors introduced by the computation. 
Discretization error is the difference between the exact solution of the PDE and the exact 
solution of the difference equations being used to approximate the PDE. The magnitude 
of the discretization error at a grid point depends on the distances between consecutive 
discrete grid points and on the number of terms in the truncated series of differences used 
to approximate the derivatives. 

Direct methods tend to be elimination methods, such as Gaussian elimination, or 
factorization methods, such as LU decomposition, QR factorization, or Singular Value De- 
composition (SVD). As seen from the systems above, most difference approximations result 
in sparse matrices, with the number of zero entries greatly outnumbering the non-zero entries. 
For this type of matrix, a method like Gaussian elimination proves inefficient, because the 
procedure fills in the zero elements with non-zero elements during the elimination, thereby 
increasing the number of non-zero elements that need to be stored in the computer and 
used in later eliminations. As a result, direct methods can be quick and accurate for certain 


matrices with special properties but prove inefficient with elliptic equations [Ref. 3]. 
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Iterative, or relaxation, methods are such that an initial approximation is used to 
calculate a second approximation which is used to calculate a third, and so on. Iterative 
methods are said to be convergent when the differences between successive approximations 
tends to zero as the number of iterations increases. Iterative methods have been much 
employed to solve Au = f because they take advantage of the sparcity of A. These relaxation 
methods are easy to implement and can be successful in determining an exact solution for 
more general linear systems than the direct methods. However, experiments have shown that 
these methods generally work very well only for the first several iterations. Convergence then 
slows and the method appears to stall. 

Iterative methods, such as the Jacobi and Gauss-Seidel methods, begin with an initial 


guess. Valuable insight can be gained by using an initial guess consisting of Fourier modes: 


oj = sin (7), 0<j<N, 1<k<N-1. (IV.12) 


The integer k is the wavenumber and indicates the number of half sine waves which make 
up the vector v on the domain. Small values of & correspond to long, smooth waves. Large 
values correspond to oscillatory waves. The rapid decrease in error during the first several 
applications of an iterative method results from efficient elimination of the oscillatory modes 
of this error. Unfortunately, the iteration is much less effective in reducing any remaining 
smooth components. Most relaxation methods possess this property, known as the smoothing 
property (Ref. 5]. The issue is whether basic iterative methods can be adapted to make them © 
effective on all error components. 

One way to improve an iterative method is to start with a good initial guess. By 


performing some preliminary iterations on a coarse grid, an improved initial guess can be 
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derived and used on the original fine grid. This coarse grid is a grid (Figure 14) with a larger 
interval between grid points. Performing these iterations, or relaxing, on a coarse grid is less 
expensive since there are fewer unknowns to be updated. 

With this is mind, recall that most relaxation schemes are effective at eliminating 
the oscillatory components of error but prove futile against the smooth components. Now, 
consider these remaining smooth components on a coarser grid. A smooth curve on a fine 
grid can appear more oscillatory on a coarse grid. Note in Figure 15 that the grid points 
on the coarse grid (bottom grid) correspond to the even-numbered grid points of the fine 
grid (top grid). Consider the kth mode of the fine grid evaluated at the even-numbered grid 
points. If 1 << k < N/2, then the components of the kth mode on the fine grid can be written 
as 


Who; = sin (2 | = sin (=) = whey. (IV.13) 


The superscripts on w indicate on which grid the vector is defined by giving the grid 
spacing on the grid. The vector w” is then defined on the fine grid, while the vector w?" is 
defined on the coarse grid. 

Thus, the kth mode on the fine grid becomes the kth mode on the coarse grid. The 
kth mode also becomes more oscillatory in passing from the fine grid to the coarse grid. 
From figure 15, k = 4 on the grid where N = 12 is 1/3 of the highest possible frequency, 
while k = 4 on the grid where N = 6 is 2/3 of the highest possible, making the wave more 
oscillatory on the coarse grid. This is true only for 1 <k < N/2. If k = N/2, the kth mode 


on the fine grid is the zero vector on the coarse grid. And for k > N/2, the kth mode on the 


fine grid becomes the (NV — k)th mode on the coarse grid. When this phenomenon, known 
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k = 4 wave on N = 12 grid 


k = 4 wave on N = 6 grid 


Figure 15. A wave with wavenumber k = 4 ona fine grid (N = 12) appears more oscillatory 
on a coarse grid (N = 6). 


as aliasing , occurs, oscillatory modes on the fine grid are misrepresented as smooth modes 
on the coarse grid. 

All of this suggests that, once convergence slows on the fine grid, the residual may 
be computed and the residual equation can be relaxed on a coarser grid to obtain an ap- 
proximation of the error. This error can then be returned to the fine grid to correct the 
original approximation. This method for improving iterative approximations, called coarse 
grid correction, leaves some issues unresolved, such as how to transfer vectors from a coarse 
grid to a fine grid (and back again), how to compute the residual on the fine grid, and how 
to compute an error estimate on the coarse grid. 

4.  Intergrid Transfers 

Before an iterative method can be used on coarser grids, a mechanism is needed for 


tranferring vectors between grids. A discussion of intergrid transfers will deal only with the 
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case in which the coarse grid has twice the grid spacing of the next finest grid, since there 
seems to be no advantage in using grid spacing with ratios other than two. 

Transferring a vector from a coarse grid to a fine grid is called interpolation, or 
prolongation. Although many methods exist, for most multigrid purposes linear interpolation 
is quite effective. The linear interpolation operator will be denoted as J},. The operator 
takes coarse grid vectors and produces fine grid vectors. With v?" as the coarse grid vector 


and v" as the fine grid vector, the notation J, v7" = v" transfers the vectors such that 
v3, = us oN 

har = 4 (oR +0") é 

This transfer notation shows that, at avecamierel fine grid points, the values of the vector 
are transferred directly from the coarse to the fine grid. At odd-numbered fine grid points, 
the value of the grid point on the fine grid is the average of the adjacent coarse grid values. 


For a two-dimensional problem, the interpolation operator is defined in a similar 


manner. According to the notation J4,v?" = v", the components of v” are given by 


h = 2h 
U23,2) = Uj 
h = 1(,,2h 2h 
V9i41,23 = a (Uji 3 + UfF) _ 
0 < 4,7] < 2 —l 
h = 1(,,2h 2h 
Voi2j41 = a(Upia1 + Uf5) 
h _ 1 (,,2h 2h 2h 2h 
Upsiaje = gUeige + lig + Uige + Uy) 


Transferring a vector from a fine to a coarse grid is called restriction. One restriction 


operator, injection, is defined by [?*v* = v2", where the components of v”" are given by 
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With injection, the coarse grid vector takes its value directly from the corresponding fine 
grid point. 
Another restriction operator, known as full weighting, is defined by ey 


where the values of the coarse grid vector are a weighted average of values at neighboring 


fine grid points, or 


= (IV.14) 


w/z 


1 
2h _ h h h : 
U; = 4 Cy + 29; | is) ’ I < J < 


The two-dimensional full weighting operator is defined similarly, where the values of the 
coarse grid vector take on a weighted average of values at neighboring fine grid points. 


Thus, 
ae ae, h h h 
Vij = T6 Caer + V9;_1,0541 1 Vaig1,2j—1 1 V2i+1,25+41 


h h h h 
+2 (vi ost + U3; a341 + V9j-1,25 F See 
i oe oo IN 
+4vh, »,| ) l < 2,7 < 2 al l. 
5, Coarse Grid Correction 


A coarse grid correction scheme may now be introduced. Coarse-grid correction is 


described as a two-grid process as follows: 


-eRelax v; times on A*u* = f on 2" with an initial approximation vi". 


eCompute the residual on the fine grid, r* = f* — Av", transfer the residual from the 
fine grid to the coarse grid, r?* = [?"r", and solve the residual equation, Arhe2h — yeh 
on the coarse grid 2?" with an initial guess e** = 0. 


eCorrect the approximation obtained on 2" (v" in the first step) with the error estimate 
obtained on 27": v* + v* + I}, e7". 


eRelax v> times on A*u* = f on 2" with the improved initial approximation v". 
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The integers v, and v2 are parameters in the correction scheme which control the number of 
relaxation sweeps before and after the coarse grid correction. Experiments have shown that 
effective choices for y; and v2 are often 1, 2, or 3. 

A*" is the coarse grid representation of the original matrix A. The original problem 
could be directly discretized on 2?" using the same techniques used to generate A”, the fine 
grid operator, but at a spacing of 2h, resulting in A?". An alternate approach is to note that 
if the error on the fine grid, e* = u” — v", lies entirely within the range of interpolation, 
denoted by R(J%,), then the error can be represented by e* = J, u?", for some coarse grid 


vector u2". The residual equation on 92" can then be written as 
Nea air =r". 
Applying Seen operator [?" to the equation 
Ahh yeh — yh 


gives 
[2h Ahh y2h — p2hyh 
which leads to the natural definition 
A ea A. 
The coarse grid equation becomes 
Per) eee 


This argument is based on the assumption that the error e” is entirely within the range of 


interpolation. This may not always be the case. If the error were entirely within this range 
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of interpolation, then solving the residual equation on the coarse grid would give an exact 
solution. Regardless, this approach gives a definition for the coarse grid operator A". This 
definition is known as the Galerkin condition. This coarse grid operator A*” may not give the 
same A”* as discretizing at twice the mesh spacing. However, the Galerkin condition may 
be used to prove such things as convergence theorems, making this condition an effective 
theoretical tool. 

In the above coarse grid scheme, there is an efficient way of solving the coarse grid 
problem Ae?" = r?4, The coarse grid problem is not much different from the original 
problem, and can be solved by the same coarse grid correction scheme. Applying the coarse 
grid correction scheme to the residual equation on the coarse grid requires a move to the 
next coarsest grid, 24", for the correction step. This process is repeated on successively 
coarser grids until a direct solution to the residual equation is economical. The resulting 
algorithm is applied recursively down to the coarsest grid, sometimes a single interior grid 
point. Corrections are then applied successively through finer and finer grids, returning 
eventually to the finest grid. This algorithm is referred to as a V-cycle , named for the 
pictoral representation of the schedule of grid visits (see Figure 16). It is defined by vig 


MV*(v*, £*) and is summarized here using a compact recursive definition: 


eRelax v; times on A*u* = f on 9 with an initial approximation v". 


elf * is the coarsest grid, then proceed to the last step. Otherwise, compute the 
residual on the fine grid, r* = f* — Av", transfer the residual from the fine grid to 
the next coarsest grid, r2* = J?*r*, and set the initial guess for 2?" to zero, v*" < 0. 
Return to the first step, v2"? — MV?"(v2", r2"). 


eCorrect v® < v" + Ii,v?" 


eRelax v> times on the original equation, A*u* = f*, with an initial guess v’. 
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Figure 16. Schedule of grids for the V-cycle and the FMV scheme. 


Recall that one way to improve iterative methods is to start with a good initial guess. 
The V-cycle uses only coarse grid correction in solving the residual equation and does not 
make use of any method to improve the initial guess. The arched of nested iteration uses 
coarse grids to obtain improved initial guesses for fine grid problems. In order to obtain 
an improved initial guess for the first fine grid relaxation of the V-cycle, nested iteration 
suggests solving this problem on the next coarsest grid. To obtain an improved initial guess 
on this coarse grid, say 127", it is necessary to solve the problem on 2*". This pattern results 
in a recursive path that leads to the coarsest grid. 

The full multigrid V-cycle (FMV) joins nested iteration with the V-cycle, resulting in 
a simple but powerful algorithm. Each V-cycle is preceded by a smaller V-cycle that provides 
the best initial guess possible. The full multigrid V-cycle, defined by v"’ + FMV*(v*, f*), 


is Summarized here: 
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elnitialize f", r2", r#*, ...; v®, e2* 


.. to zero. 


eRelax on the residual equation on the coarsest grid (or solve directly if the coarsest 
grid is small). 


eCorrect e#* + e* + [4he® and perform a V-cycle using the grid 04" as the finest 
grid and 92°" as the coarsest grid. This is abbreviated as e** + MV**(e*, r*). 


eCorrect e?* using the result from the previous V-cycle, e** <- e+ I Ze 

eet! MV 2h(e2h, 2h), 

ev e vis TR, 

ev’ & MV*(v", f*). 
The last step .in the above procedure is the original V-cycle discussed earlier, but is preceded 
by several smaller V-cycles (Figure 16) in order to make an improved initial guess. 

The full multigrid V-cycle brings together many ideas and techniques, that, when 

taken alone, can have serious limitations. Full multigrid integrates these ideas so that they 


may work together to remove these limitations. 


B. THE FAST ADAPTIVE COMPOSITE GRID METHOD 

The fast adaptive composite grid method (FAC) is a natural extension of conven- 
tional multigrid methods applied to problems discretized on composite grids. Multigrid was 
developed as a global grid solver but is not efficient for cases where fine grid resolution and 
accuracy are needed only in a local subdomain. While a certain amount of accuracy may be 
suitable to correctly solve the physical system in most of the global domain, greater accu- 
racy is most often needed in a given local region. Obtaining the desired accuracy may take 


a resolution with a mesh size A, which is smaller than that in the rest of the global domain. 
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If a global fine grid were incorporated, unnecessary computations would result outside the 


local region, often at a great computational cost [Ref. 4]. 






Figure 17. A global domain Q. QQ» 1s the local domain with a smaller grid spacing while QQ; 
is that portion of the global domain Q) that retains the original coarse grid spacing. 


To arrive at FAC, the desired method, ideas will be drawn from several simpler, less 
efficient methods. It should be noted here that this development will follow closely with that 
of reference [Ref. 4]. 

1. Local Multigrid 

A natural way to avoid using a global fine grid without sacrificing the necessary fine 
grid resolution in the local subdomain is to modify the standard multigrid scheme by simply 
suppressing relaxation outside of the subdomain. Let Gb, represent a restricted relaxation 
process on the local subdomain , (Figure 17). Let the decomposition of u” be represented © 


by 
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where u® represents the unknown values of u found on the local subdomain, $b, while 
the remaining unknown values, u”, are found on the complement of the subdomain, Cao. 


Applying a similar decomposition to f”, the restricted relaxation operator has the form 


h 
h (.h gh u 
Go, (ug, fy’) a ( Gut, fh) 
where G" is a relaxation scheme corresponding to the local problem on (22, resulting in 


changes to u® only. The local multigrid method , defined by u" <- MG6, (u"*, f"), can then 
be expressed in the two-grid form as follows: 

of? g J2h(gh — Ahyh) 

eou2*  (A2h)-1 42 

eu' cu? + Ju 

eu® & Gb, (u*, f*) 
The objective of the local multigrid scheme M G%,, is to eliminate unnecessary computations 
on CQ. But only the relaxation process has been suppressed there, not the intergrid trans- 
fers or the residual calculations. Fortunately, the intergrid transfers and residual calculations 
can be suppressed without changing the final result. With u” as the initial approximation 
for one MG, cycle, u*" as the correction computed in the above two-grid form, and vi 
as the change in u* + J/,u?" resulting from the relaxation on 3 (again from above), the 
approximation produced by MG%, is u’+J},u**+v". Since va is zero on CQ2 (no relaxation 
occurs on CQ, in the relaxation scheme), J#,A,v" is zero on C4, or rather on the set of 


points that are not in 92 or its interface ( N#(|CQ2). Using the two-grid form above and 
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the Galerkin condition [?* A’ J?, = A" results in 





hgh — AM? 4 Tu??4v")) = TpP(f* — Aha") — IRMA" u?* — Th Ab" 
_ fag (f" ” A*u’) — A2hy2h _ Pe Abyh 
= frh - f2h a [Ph Ay" 


= a1 Av". 


This implies that the transferred residual f?" in the next cycle of MG%, will be zero at the 
points of CQ. In general, local fine grid relaxation has no immediate effect on the coarse 
grid equations at the points on CQ}. There is no need to use the fine grid in that region 
[Ref. 4]. 

2.  Bordered Multilevel Method 

To prevent unnecessary work away from the local fine grid, the role of u2* in CN} 
must change from representing error corrections to representing the current approximation. 
A modification of the local multigrid scheme MG%,, uses this u2" to accumulate corrections. 
Also, when using this modification of M Go, the fine grid must have enough data to compute 
1?h(f" — A*u') at the coarse grid interface points. This requires interpolation on 26 and on 
the two grid lines bordering (22. Figure 18 shows the first border and the extended local fine 
grid, 0%, , which includes both borders. 

These ideas provide for a new algorithm, called the bordered mulilevel scheme, that 
avoids all fine grid work away from the local fine grid 02. In this new algorithm, residuals are 
computed on the extended local fine grid, which includes the local fine grid, the boundary 


points of the local fine grid, and the first border. The residuals are transferred to the 
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& first border point 


@ local fine grid point 


Figure 18. The local grid Q2, its first border (indicated by ¢), and the extended local grid 
Oh 


coarse grid points lying under the fine grid and the interface (residuals are also computed 
on the coarse grid outside the fine grid). A global coarse-grid correction is computed and 
interpolated to the extended local fine grid. Relaxation is then performed at fine grid interior 


points. 


3. Multilevel Composite Grid Method 

The bordered multilevel scheme requires two additional grid lines bordering the local 
fine grid in order to compute the residual correction at the coarse grid interface points. This 
residual correction is the composite grid residual at the coarse grid interface points. Let 
2 denote the composite grid, M2 = N?*U". By generating a composite grid operator, 
A& : U4 — U4 where U® is the nodal space associated with the composite grid M4, the 


residual correction can be computed more directly without the need for any borders. This 
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leads to the multilevel composite grid method , represented by u# + ML4(u, f°), and defined 
as follows: 


eSet £2" = 17°(f4 — A®u*). 


eCompute u2” < (A2")—! £2", 


eCompute uf + ué + [3 u*. 


eSet f° = Ip(f4 — A’u®). 
eSet u* — 0 (the initial guess for the fine grid is zero) and compute u” + G"(u", f"). 
eCompute ut ~ ué + JAu’. 

In developing this multilevel composite scheme, it has become necessary to define 


some new notation. First, the domain {2 needs to be partitioned as 
2 =X Ur Qr (IV.15) 


which breaks the domain into subdomains. Qc is that part of the domain that contains all 
coarse grid points excluding the interface points, (22; includes the interface points, and !20r 
includes the fine grid points. On the composite grid in Figure 19, the grid 24 is partitioned 
into the following: 

Om = 22 JOFU OG. (IV.16) 
ne consists of the coarse grid points (excluding the interface points), Gi includes the coarse 
grid interface points, and O52 is the fine grid patch. 

A block decomposition of the operators and grid equations results from the partition 


in (IV.16). The block decomposition of the grid equations A*u4 = f2 is given by 


h h hk h he 
Acc Ac: Acr Uc fc 
h A h h _ A 
Aico Ay Alr I = f; 
h h h h 0% 
Apc Ap, Apr Up fr 
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© Coarse grid point 
LJ Interface point 


é Fine grid point 





Figure 19. Composite grid partition. Subgrid ne contains only coarse grid points, subgrid 
Qt contains only interface points, and subgrid O. contains only fine grid points. 


The notation AS signifies that the equation being solved corresponds to an X-grid 
point and the stencil reaches to a Y-grid point. If 4, and J?" represent the grid transfer 
operators restricted to the local fine grid, then certain observations result from the way the 
grid operators are constructed. For the transfer equation ut = J4 u?*, the transfer operator 


is 


10 0 
m = | O70 |, (IV.17) 
Ox Ih 


which may be interpreted as follows: outside the local fine patch, interpolation is simply 
the identity. Interpolating u** inside the local fine patch is accomplished using standard 
prolongation away from the borders. The star, x, represents actual storage of the slave 


points as ne values. 


of 





The transfer operator J?” is given by 


10 0 
Lo =. 0 Eb & (IV.18) 
00 IP 


Residuals outside the local fine patch are transferred by injection to correct the equations on 
the coarse grid (02? U2"). These residuals do not affect Nf equations. Rather, residuals in 


ne he ee 
= are used to correct 22" equations. Residuals in 27 UM are used to correct 27" equations. 
F F &q I F I 


The transfer operator J; is given by 


0 
Ip = | 0 |. (IV.19) 
I 


Residuals are transferred from the composite grid to the fine grid patch in the trivial manner. 


The transfer operator ibs is given by 


if = (0 0. 7): (IV.20) 


Interpolation from the fine grid patch to the composite grid is the trivial interpolation, the 


identity. 


The composite grid operator is given, as before, by 


A h h 
cc Gl “cr 
Ae = 1C Ary Air (IV.21) 


bh bh & 
Apo Ap, Apr 
Recall that the notation 42, means that the equation being solved corresponds to an X-grid 


point and the stencil reaching to a Y-grid point. Since the coarse grid stencils do not reach 


. h h 
any fine grid points, Apo and AgGp are zero. 
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The composite grid operator can then be written as 


Ar, x 0 
AB = | .® x & |. | (IV.22) 
0 -* Af 


The composite grid operator is block tridiagonal and agrees with the coarse grid operator on 
' Qc and with the fine grid operator on Qr. The stars, x, are non-zero but unspecified entries 
whose values are not important at this time. 

The composite grid residuals in Qc and Qf can be computed using the respective 
coarse and fine grid stencils. Special evaluation is, however, required at the coarse grid 
interface points. The stencils for A* in the multilevel composite grid scheme M L* must be 
computed for use on Oe. 

The multilevel composite grid scheme ML given above is written in the immediate 
correction form. This scheme is considered to be “immediate correction” because both 

h 


intermediate quantities u2" and u” attempt to approximate the current algebraic error e* = 


u2” — u4 (u”’ is the exact solution) while solving the residual equation 
Afe® = f2— A*ut*. (IV.23) 


This immediate correction form attempts to solve the residual equation (IV.23) by trans- 
ferring the composite grid residual to each level in turn, computing a correction there, and 


immediately interpolating it back to the composite grid. 


4. The Fast Adaptive Composite Method 
A final modification to the multilevel composite grid scheme M L/ is made by replacing 


fine grid relaxation with a direct solver. The replacement of the fine grid relaxation with a 
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direct solver is given by 
G"(u", f") = Chg aes in 
This modification gives the fast adaptive composite grid method (FAC). This method is 


represented by u2 + FAC#(u4, f2) and is defined as follows: 


eSet £7" = 7*(f4 — ASut). 
eCompute u2" + (A2")—1 £2", 
eCompute u4 + ub + 2 u, 
eSet f* = Ih(f — Abut), 
eCompute u” «+ (L")—!£". 


eCompute u2 + uf + [tu’. 


C. DISCRETIZATION ON THE COMPOSITE GRID 

To apply FAC, the equations to be solved must be discretized on a composite grid. 
FVE provides the method for this discretization since FVE has proven effective with PDEs 
on a composite grid. 


1. Discretization of the Model Problem on a Composite 
Grid , 


As an example, recall that the model problem is the potential equation 


V-(pVv) = n in (IV.24) 
y = W on ODy Ug (Dirichlet) 
(pVwv)-nm = yy on OQwUONs (Neumann). 


The composite grid can be taken as a coarse grid with a patch of fine grid points, where the 


coarse grid points and the slave interface points along the boundary of the fine grid patch 
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appear in the patch equations as if they were Dirichlet boundary points for the patch [Ref. 
4]. The patch volumes are chosen first, with the volumes along the interface being largely 
dictated by the neighboring patch volumes. Triangulation follows with the patch geometry 
being determined at the interface by combining pairs of triangles of opposing orientation 


that abut the slave points (Figure 20). 





Figure 21. Interface point volume for a composite grid. 
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As an example of a stencil for a composite grid, consider the interface point P in 
Figure 21. The control volume is smaller for the interface point than a control volume for 
the same point on a uniform coarse grid, and there are now seven triangular regions over 
which to integrate. In the seven regions are ten portions of the border of the volume. Using 
the same procedures as for the general interior grid point of the uniform grid, integrate over 
each of the ten segments of the volume boundary. The resulting stencil form of the equation 
for the interface point P is 


Pi5+4 


Here, as before, the notation © represents the sum of the surrounding stencil entries. 


2. Discretization of the Weldpool Problem on a Composite 
Grid 


With an understanding of how F'VE can be used to discretize a PDE on a composite 
grid, now consider discretization of the weldpool problem on a composite grid. Equations 
II.5 - II.8 are considered now on a composite grid. Figure 22 shows a two-level single patch 
composite grid, with grid spacing on the fine patch taken as half the grid spacing of the 
coarse grid. 

The composite grid points that will be examined are those along the coarse-fine inter- 
face. Six grid points are labeled in Figure 22. These grid points are interior interface points 
(P and Q), Neumann boundary interface points (M and N) and Dirichlet boundary interface 


points (D and F). For the Dirichlet boundary, recall that since the Dirichlet conditions are 
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Figure 22. A simple composite grid with volumes around interior boundary points P and 
Q, Neumann boundary points M and N, and Dirichlet boundary points D and F (extended 


volumes around D1 and F1). 


imposed directly on T, the control volumes are extended so that the control volumes nor- 
mally associated with the boundary points are taken into account by an extended volume 


from the neighboring points (D1 and F1). 


Beginning with the conduction equation for the solid, 


= = Ma VT, — (IV.25) 


the discretization of the conduction problem resulted in the following integrals: 


ue 1 n — n a 
S- F I gey'dV — aMa , [ Vers! ads] 1s (IV.26) 


i,j=0 


= = 1 n a —} now n 
= > = [av +0 o)Ma [ ver, ads T™. 


ij=0 
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As before, the integrals over the basis functions are calculated by representing the basis 
functions as collections of triangular planes. The triangular planes are assembled to form 
the “hat” function. An interior grid point along the interface has seven triangular planes 
joined in this manner (Figure 21). A Neumann boundary point has four such triangular 
planes, as can be seen in Figure 23. A neighboring Dirichlet boundary point, with its 
extended volume, also has seven triangular planes. (Notice that the edge labeled “SZ3” 


does not lie under the “hat” function of the neighboring Dirichlet boundary point D1.) 





Figure 23. Neumann boundary point (M) volume and Dirichlet boundary point (D) extended 
volume on a composite grid. 


Following the techniques outlined above for a uniform grid, interface point stencils 
have been computed for the volume and surface integrals of equation (IV.26). Referring to 
Figure 22, the stencils have been computed for interior interface grid points (P and Q) on the 


top and bottom boundary, and also for the Neumann (M and N) and neighboring Dirichlet 
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(D1 and F1) grid interface points on the top and bottom boundary. Below are the stencils 


for the volume and surface integrals for point P: 


| 32 — 5e 16 + be 
h?r 
[ (SX Thadea] dV = 7 | 32 1e 192 + 5e 16 + 6e | Ti; 
Yi Ug 32—10e 48+4€ 16+ 6¢ 
and 
32 
| (5 Tra V Prt 2 dS = = 16 — 9e —128 16+ 9e | Tj. 
S35 kg 16—6e 32 16+6¢ 


In the above stencils, « = A (r is the radial coordinate for the central point of the 
volume), and, as before, the 27 resulting from the integration in the @ direction has been 
factored out. The stencils for the remaining boundary points can be found in the Appendix. 

The next stencils to be calculated are those belonging to the convection-diffusion 


equation. Again, observe that the discretization of the convection-diffusion equation, 


ea +V-(uT) =Ma'V°T, (IV.27) 


has the same form as the conduction equation for the solid except for the term 
I V - (uT)dV. | (IV.28) 


Interface point stencils have been computed for the integral in equation (IV.28). Referring 
to Figure 22, the stencils have been computed for interior interface grid points (P and Q) 
on the top and bottom boundary, and also for the Neumann (M and N) and neighboring 


Dirichlet (D1 and F1) grid interface points on the top and bottom boundary. The stencil 
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y | 


for point P has the following form: 


where the entries of the stencil are themselves stencils applied to the function YW. The entries 


are given here: 


2 =2 
A= —4 2 Vi BS 2 Vij 
4 —2 
C' — —] Wiis D — 1 1 Wj 
—3 | —l 
2 
3 =2 
G= 1 St NW H = -2 | ¥,, 
22 —4 2 


The stencils that have been developed are for the conduction equation in the solid 
(IV.25) and the convection-diffusion equation in the liquid (IV.27). Also, only stencils for 
upper interior interface grid points are displayed here. Many stencils have been developed 


and are displayed in the Appendix. 
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3. Development of the Weldpool Problem 


The stencils developed in the previous section are for a model of the weldpool problem. 





The model used to develop the above stencils assumed a stationary boundary and considered 


only a flat interface between the two states (Figure 24). 








Figure 24. An ezample of a composite grid for the simplified weldpool problem with a sta- 
tionary, flat boundary. The grid spacing is finer along the solid-liquid interface. 





The real problem is far more complicated and requires some method for tracking the 


solid-liquid interface. The next step in solving the weldpool problem is to form the boundary 





into its naturally occuring shape and place a fine grid over the solid-liquid interface, where 
the boundary is still considered to be stationary. An example of this composite grid is found 
in Figure 25. 

The next step is to consider the movement of the interface. The movement results 


from the melting and solidifying of the solid. This movement is governed by the local heat 
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Figure 25. An example of a composite grid for the weldpool problem with a stationary, curved 
boundary. The grid spacing is finer along the solid-liquid interface. 


balance around each interface point. A requirement for a control volume along this interface 
is that both the current time and the next time step must be contained in the volume. Not 
only the current position of the interface, but an approximation of the interface after one 
time step are required to construct the current local grid [Ref. 1]. This will require having 
an adaptive component in the time stepping. That is, at each step, the next position of the 
interface must be computed. If the change in position is too large and the new interface 
point lies outside the control volume, then the time step must be made smaller. A new local 
grid must also be computed at each time step. In this way, the control volumes are chosen 
in a manner that allows the interface to be moved. 

These, and other complicating factors, indicate that a good deal of work remains 


before the full weldpool problem can be considered “solved”. 
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Vv. CONCLUSION 


The research for which this work is a part is designed to develop efficient 
numerical methods for solving the full weldpool problem. This work was undertaken with 
the goals of developing the FVE method for spatial discretization and finite differences for 
time differentiation to discretize the equations that describe the weldpool problem, and 
develop stencil equations for use in FAC. This task began with the application of finite 
differences to the time derivatives in the conduction problem, resulting in a time-stepping 
scheme that makes use of a weighted average at the current and the next time steps. F'VE is 
then applied to the potential flow problem to develop the method for discretizing the spatial 
derivatives. Having developed these methods on model problems, the weldpool problem is 
tackled by applying finite differences and FAC to discretize the governing equations. 

Multigrid is presented, including a development of the coarse grid, intergrid transfers 
and coarse grid correction. Multigrid is introduced here because of its effectiveness as an 
iterative method. Multigrid proves sufficient for solving problems on a uniform grid. How- 
ever, the weldpool problem requires an extremely fine level of accuracy along the liquid-solid 
interface. To achieve this level of accuracy on a global grid would be too costly. Thus, FAC 
is developed to make use of composite grids. These grids allow finer grid spacing along the 
boundary of the molten pool and the solid, while still providing for coarser grid spacing 
where fine accuracy is not required. To arrive at FAC, several ideas are drawn from sim- 
pler, less efficient methods related to multigrid, including local multigrid and the multilevel 


composite grid method. A simple model of the weldpool problem is then discretized on a 
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composite grid using FVE. Stencil equations are developed for the conduction equation in 
the solid and the convection-diffusion equation in the liquid. 

Finite Volume Element methods are shown to be an effective method for discretizing 
equations on a composite grid. The stencil equations that result from the discretization will 
be coded into FAC at a later date to determine their effectiveness. Future work will center on 


applying the FVE discretization and FAC to a moving, time-dependent liquid-solid interface. 
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APPENDIX. STENCIL EQUATIONS 


Included here are the remaining stencil equations that were not included in the main 
body of the text. 
The first of the stencil equations are for the conduction problem. The discretization 


of the conduction problem resulted in the following integrals: 


3 F . oldv -— aMa” [ Vert! - nds] Tr (1) 
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Figure 26. A composite grid with volumes around interior boundary points P and Q, Neu- 
mann boundary points M and N, and Dirichlet boundary points D and F (ertended volumes 
around D1 and F1). 
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Referring to Figure 26, the stencils have been computed for the conduction problem 
above (.1) at interior interface grid points (P and Q) on the top and bottom boundary, and 
also for the Neumann (M and N) and neighboring Dirichlet (D1 and F1) grid interface points 


on the top and bottom boundary. Given here are the stencils for the volume and surface 


integrals for points Q, Dl, F1, M, and N. 


At point Q: 
pee 16—6e 48—4e 32+ 10¢€ 
/ ( Tb dV: = 384 16 — 6¢ 192 — de 32+ lle Lia, 
ti \ kl 
16 — de 32 + OE 
At point D1: 
32 — O€ 
h?r | 
J (Studs) av = aay | 32 Llc 224 + 26¢ i 
"7 \ kd 32—10e 48+4¢€ 48+ 30¢e 
At point FI: 
ee l6—6e 48—4e 48+ 20€ 
i ( Tutus dV = 384 16 — 6€ 208 + 5e dis: 
‘i \ kl 
16 — 5e 48 + 49e 
At point M: 
S+eE 16 + 5e 
h?r 
| Tides) Vo = 104 + 24¢ 16 +6¢ | Ti. 
Yi kg 32+6e 16+6¢ 
At point N: 
her 16+2€ 32+10¢ 
I (5 Ton dV = 384 88 + 19€ 32 + lle T; % 
Wi VE 


24 + 6e 
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At point Q: 


= —| 16-9 ~128 16+9¢ | T;;. 


16-—6e 32 16+6¢ 
dS 
32 


/ (5 TV r+ 2 
Sui ku : 
32 


At point D1 (Notice here that T;,; is known for the furthest right entries of this stencil. 


These three stencil elements can be pulled from the stencil.): 


32 16 + 12 
[ (5 TV PRL s dS = = 16 — Ye —128 — 9e —32 — 18€ | T;,;. 
Md ded 16 — 6e 32 16 32+ 306¢ 


At point F1 (Notice here that 7;,; is known for the furthest right entries of this stencil. 


These three stencil elements can be pulled from the stencil.): 


‘ 16 — 6€ 32 16 32+ 30¢e 
I (5 Th V Ons: | dS = — | 16—9e —112 + 7e —16 —10€ | J;,;. 
Sij k,l 32 | 
16 — 6€ de 
At point M: 
16 + 4e 
/ (5: TkV Prt 2 dS = 7 —64 — 2le 16+ 9e | T;;. 
ms \ ke 16+2e 16+6¢ 
At point N: 
; 164+ 2e€ 16+ 66¢ 
| (5: TeV Pr ° i dS = = —64 — 2le 164+ 9e | T;,;. 
Si,j kl 32 
16 + 4e 


In the above stencils, € = 2 (r is the radial coordinate for the central point of the 
volume) and, as before, the 27 resulting from the integration in the @ direction has been 


factored out. 
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Discretization of the convection-diffusion equation, 


a + V+ (uT) = Ma 'V’°T, (.2) 


has the same form as the conduction equation for the solid except for the term 


[ V -(uT)dV. (.3) 


Interface point stencils have been computed for the integral in equation (.3). Again referring 
to Figure 26, the stencils have been computed for interior interface grid points (P and Q) 
on the top and bottom boundary, and also for the Neumann (M and N) and neighboring 
Dirichlet (D1 and F1) grid interface points on the top and bottom boundary. The stencils 
given here ard for points Q, D1, Fl, M, and N. 


At point Q, the stencil has the following form: 


ABC 
[v-(wtav -~—|D EB FIT, 
V 
H 


G 


where the entries of the stencil matrix are themselves stencils, applied to the function W. 


The stencil entries are given as: 


2 —4 2 2 
A=| -2 Wi B=| -1 1 Wi; 
—2 3 
C’ = —] 3 Wes p= —] Wi; 
—2 
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—3 
EB = l 1 Wi FP — —] Wij 
9 4 
G=|2 Wii H= 2 —-4 |V;, 
—2 2 
At point D1, the stencil has the following form: 
A 
/ Yu se — |B -e f. 
: SS a 449) 
V 16 DEF 


where the entries of the stencil matrix are themselves stencils, applied to the function W. 


The stencil entries are given as: 


2 4 
A=|-4 2 Wi B= ae Wi 
—3 
—2 4 

C=|1 Wi; D=|3 <1 va 

1 —4 —2 
E= i i F = 3 -3 |W, 

22 -—4 2 —2 
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where the entries of the stencil matrix are themselves stencils, applied to the function WV. 


The stencil entries are given as: 


Z ae ae 
A=] ~2 Vi, B=| -1 1 Wi 
me, oe 3 
C' — —3 3 Wi D= —] Wi 
—2 
—] 4 
= 1 Wi = 2 Wi 
—4 =2 
G= 4 4 ]W,, 
Z —2 


At point M, the stencil has the following form: 

A B 
C D kee 

EOF 

where the entries of the stencil matrix are themselves stencils, applied to the function W. 


The stencil entries are given as: 


a —2 4 
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—4 | 2 


At point N, the stencil has the following form: 


A B 
[/v-(nav = = C DIT, 


E 


where the entries of the stencil matrix are themselves stencils, applied to the function W. 


The stencil entries are given as: 


2 9 —2 
A= Wij B= —] 3 Wij 
eA. <j —3 
C= 1 |, D= ~3 21, 
2 4 
E= 2 —4 Wij 
2 


The stencils that have been included here are for the conduction equation in the solid 


(IV.25) and the convection-diffusion equation in the liquid (IV.27). The stencils here are for 
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the interface, either in the interior or along the boundary. Stencil equations for the vorticity 


and stream functions, equations (II.7) and (II.8), still remains for future computation. 
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